#include "Multigridsolver copy 2.h"
long double basic_f2rd(long double * P)
{
    return -(sin(P[0])-cos(P[0])*cos(P[0]))*exp(sin(P[0]));
};

long double basic_f(long double *P)
{
    return exp(sin(P[0]));
};
int main()
{
int n=128;
Function _f=basic_f2rd;
long double _a[1]={0};
Multigrid_1D_Solver case1;
case1.Initialize(n,5,_f,1,exp(sin(1.0)));
case1.setRestircType(1);
case1.setIntpType(1);
case1.setSmoothType(2);

VectorXd u_,f_,r_,r_origin,acu;
acu.setZero(n+1);
for(int i=0;i<=n;i++)
{   
    _a[0]=1.0*i/n;
    acu(i)=basic_f(_a);
};
long double acu_norm=acu.lpNorm<Eigen::Infinity>();
long double h;

u_=case1.get_u_origin();
f_=case1.get_f_origin();
h = case1.get_h();
r_origin=case1.ResidualSolver(u_,f_,h);
double r_origin_norm=r_origin.lpNorm<Eigen::Infinity>();
int cnt;
for (cnt=1;cnt<=500;cnt++)
{
    u_=case1.Vcycle(u_,f_,h);
    r_=case1.ResidualSolver(u_,f_,h);
    if (r_.lpNorm<Eigen::Infinity>() <= r_origin_norm*2.2e-16)
    break;
};

// int cnt;
// for (cnt=1;cnt<=1;cnt++)
// {
//     u_=case1.FMGcycle(f_,h);
//     r_=case1.ResidualSolver(u_,f_,h);
//     wucha=u_-acu;
//     if (wucha.norm()/acu_norm <= 1e-3)
//     break;
// };

//  std::cout<<"The residual is\n"<<r_<<std::endl;
//  std::cout<<u_<<std::endl;
std::cout<<"Residual="<<r_<<std::endl;
 std::cout<<"cnt= "<<cnt<<std::endl;

std::cout<<"Residual norm_inf="<<r_.lpNorm<Eigen::Infinity>()<<std::endl;

std::cout<<"relative accu= "<<r_.lpNorm<Eigen::Infinity>() /r_origin_norm<<std::endl;
// for(int i=1;i<2;i++){
//    case1.Vcycle(); 
// }
// solution1=case1.Vcycle();
// // solution2=case2.Jacobi();
// // acu.resize(solution1.rows());
// for(int i=1;i<n;i++)
// {
//     // acu(i-1)=exp(sin(1.0*i/n));
//      acu(i-1)=sin(M_PI*i/n);
// };
// long double p=(acu).norm();
// std::cout<<solution1<<std::endl;
// std::cout<<acu<<std::endl;

// // std::cout<<"Jacobi:"<<std::endl;
// std::cout<<((solution1-acu).norm())/p<<std::endl;

};